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Abstract 

Transport models of growth hormones can be used to reproduce the 
hormone accumulations that occur in plant organs. Mostly, these accu- 
mulation patterns are calculated using time step methods, even though 
only the resulting steady state patterns of the model are of interest. We 
examine the steady state solutions of the hormone transport model of 
[Smith et al| p006| for a one-dimensional row of plant cells. We search for 
the steady state solutions as a function of three of the model parameters 
by using numerical continuation methods and bifurcation analysis. These 
methods are more adequate for solving steady state problems than time 
step methods. We discuss a trivial solution where the concentrations of 
hormones are equal in all cells and examine its stability region. We iden- 
tify two generic bifurcation scenarios through which the trivial solution 
loses its stability. The trivial solution becomes either a steady state pat- 
tern with regular spaced peaks or a pattern where the concentration is 
periodic in time. 

keywords: Bifurcation analysis pattern formation parameter depen- 
dence auxin transport model stability periodic solution pattern 
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1 Introduction 



1.1 Biological background 

For centuries, the formation of well-defined patterns in plants, such as the ori- 
entation and shape of leaves, their venation patterns, the spatial distribution of 
hairs and stomata, the early embryonic development patterns and the branch- 
ing patterns in both root systems and treetops, has intrigued many scientists. 
Experimental research has identified a number of molecular components that 
play a major role in several of these pattern formation processes. One of them 
is the plant hormone auxin, and more specifically the auxin molecule Indol-3- 
Acetic Acid (lAA). Experiments have shown that the active directional trans- 
port, which leads to accumulation spots of the auxin hormone, plays a central 



borough et al 2011) 



part in the pattern formation (Scarpella et al 2006 Benkova et al 2003 Bils 



Based on such experimental evidence, Reinhardt et al. developed a con- 



ceptual model that describes the auxin transport through the cells (Reinhardt 



et al 2003). Smith and collaborators then constructed a computational simula- 



tion model ( Smith et al 2006 ) incorporating the experimental evidence that the 
transport of the auxin molecule lAA is driven by a pumping mechanism that is 
mediated by PINl proteins located at the cell membrane in addition to diffu- 



sion ( ^Palme and Galweilcr^ 1999). Therefore, Smith and collaborators modeled 
the transport of the lAA hormone through the cells by describing the simulta- 
neous evolution of the PINl protein and the lAA hormone concentrations over 
time. Also other computational models were developed based on these molecular 
mechanisms identified by Reinhardt et al. For instance Jonsson et al. proposed 



a phyllotaxis model based on the polarized auxin transport (Jonsson et al 2006 ) 



They analyzed a simplified version of their model that assumes an equal and 
constant PINl concentration in every cell and membrane. In their simulations 
they used a linear row of uniform cells with periodic boundary conditions. The 
results show that the spacing and the number of peaks in this simplified model 
depends on the different parameters. Jonsson and collaborators also performed 
a stability analysis and found an analytical expression for the eigenvalues, be- 
longing to a solution pattern with equal auxin concentrations. The eigenvalues 
are all real and a function of the model parameters of the model. They also iden- 
tified the parameter threshold where the largest eigenvalue becomes unstable. 
Beyond this threshold, all solutions will contain auxin peaks. 

This paper expands the study of the steady states in the transport of hor- 
mones. We limit ourselves to the study of the auxin distribution in a linear row 
of uniform cells that represents, for example, a cross section through a young 
leaf. We perform a thorough mathematical exploration of the behavior of the 
models and how their equations are solved starting from the coupled model of 
Smith et al (2006). In contrast to the analysis of Jonsson et al. on a row of 
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cells where the concentration of PINl was fixed, we will use a coupled model 
where the PINl concentration is allowed to change from cell to cell. The anal- 
ysis gives new insights into the spacing of auxin accumulations that form the 



basis of vascular development ( Scarpella et al 2006 1 



The main contribution of the paper is a systematic numerical bifurcation 
analysis for the coupled model describing the transport of growth hormones. 
The analysis identifies two generic bifurcation scenarios that reappear for vari- 
ous choices of the parameters of the problem. Through the bifurcation diagrams 
we identify the genesis of the patterns that were observed by Smith and col- 
laborators. Furthermore, we have found a limited parameter range that allows 
periodic solutions in the system. In these solutions, the concentration of auxin 
of each cell varies periodically over time. To the knowledge of the authors these 
results have not appeared in the literature. 

We present our work as follows. In section |1.2| the basic cell polarization 
and auxin transport model of Smith et al., where PINl is allowed to change, is 
reconstructed using graph theory to describe the geometry and connectivity of 
the cells. In the next section, a specific model that will be used in the simulations 
is defined. In Smith et al ( 2006 ) the model roughly correlates to a ring of cells 
around an axial plant organ. In subsequent situations we consider a linear row 
of cells running from the margin to the midvein of the leaf and we consider 
the fiuxes at the boundary of the leaf. We will describe this by using Neumann 
boundary conditions instead of periodic boundary conditions which is explained 



in section 2.1 Also in this section we look at the different parameter values. 
Similar to Smith et al, we use time integration to solve the coupled model of 
Smith et al. in section |2.2| Since we are only interested in the steady state 
solutions, we define in section [273} the corresponding steady state system. For 
this steady state model, we define in section [3] a trivial solution and its stability 
properties. The stability is dependent on the model parameters and we examine 
for which parameter regions the trivial solution is stable. Section [4] contains the 
techniques that will be used to solve the model. In particular, we will discuss 
bifurcation analysis (4.1 1 and continuation methods (4.2). Bifurcation analysis 
reveals the relation between the stability of a solution and the model parameters 
and continuation methods calculates approximate solutions in function of a 
model parameter. In section [5] we show the results of our simulations. In 
section |6] we conclude and give an outlook. 



1.2 Description of the mathematical model 

Before constructing a compartmental model that describes the concentration of 
growth hormones per cell and its transport through a plant organ such as for 
example a leaf, its geometry must be specified. We do this with the help of 
graph theory. 
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1.2.1 Geometry of the cells 



The two-dimensional model of a leaf consists of cells and cell walls and the 
geometry can be considered as a structure of piecewise linear cell walls. We 
represent this structure of cell walls by a graph G {E, V) with edges E and 
vertices V. In this graph the set of edges E represents the cell walls of the leaf. 
Each edge connects a pair of vertices, the endpoints of a cell wall. We label the 
vertices with the letters k and I {kj £ V) and an edge with the tuple {k,l). 
Another characteristic of a cell wall is that it has no direction. The edge {k,l) 
represents the same edge as {I, k) and we refer to it with the pair {k, I}. Because 
cell walls only intersect at their endpoints, the graph G is a planar undirected 
graph. Each cell wall has a length and therefore the graph is also a weighted 
graph. We denote these weights, the length of the cell wall {k,l}, as Iki- In the 



model of Smith et al (2006), only the length of the wall that separates two cells 
is important and not its shape. We therefore assume in this case that there 
exists at most one straight edge that separates two different faces of the graph 
G. 

A graph G that represents the geometry of a leaf is displayed with full lines 
on figure [l] A cell in this model of the leaf is a face of the planar graph G. 
We refer to a cell with the letters i and j. As a consequence a cell is a vertex 
in the weak dual undirected graph H (^E, of graph G {E, V) (displayed on 

figure [l] with dashed lines). In this graph H every vertex i £V represents a cell 
and every edge G E represents the connection between two neighboring 

cells and is therefore always undirected. Because an edge is the connection 
between neighboring cells, we can define the neighboring cells of a cell in the leaf. 
They are represented by the adjacent vertices of vertex i denoted as M C V. 
Furthermore, the dual graph H is, just like graph G, a weighted graph with 
the same weights as G. Indeed, for each edge £ E we associate one edge 

{k, 1} G E which represents their common cell wall. Thus we can associate the 
length of the cell wall Iki with We denote it with lij. 

The graphs G and H represent the geometry of a group of cells in a 2D plant 
organ at a certain moment in time. In advanced models this geometry changes 



over time because cell walls grow and cells divide ( Smith et al 2006 1 , but in this 
paper we look at the basic, coupled cell polarization and auxin transport model 
of Smith et al. and the geometry is assumed to be static. As a consequence also 
the graphs G and H are fixed. 

1.2.2 Transport of growth hormones 

With the help of the geometry it is now possible to formulate a model for the 
transport of the growth hormones through the cells. We will formulate the 



coupled basic model of Smith et al ( 2006 1 applied to a section across a leaf with 



the geometry specified above. 

In every cell i two substances play an important role in the growth process: 

• The concentration of proteins PINl in cell i, which is time dependent, is 
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Figure 1: A graph (solid line) and its weak dual graph (dashed line) that models 
the static geometry of the cells. The edges of the graph are modeling the cell 
wall and only the length plays a role in the hormone transport. These lengths 
Ihi are the weights associated with the edge between vertex k and /. Similarly, 
each edge of the weak dual gets a weight corresponding to the length of the wall 
between cell i and j. So it gets the same weight = Iki- 
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denoted as pi {t) G M+, 



• The concentration of the hormone lAA in ceU i, also known as auxin, is 
also time dependent and is denoted as (t) G M+. 

The model describes the evolution of these concentrations in each cell. This 
evolution depends in a non-linear way on the concentrations of the neighboring 
cells. The value of pi (t) is determined by the production and decay of PINl. 
Its time evolution for each cell i £ V is modeled by 

dpi{t) PpiNo +PpINa^ W 

where Ppj^^ G is the base production of PINl proteins, ppj^ G M+ is a 
coefficient capturing the up-regulation of PINl production by auxin, Kpj^j G 
is the saturation coefficient of the PINl production and /ipj„ G M"*" is the PINl 
decay constant. This means that the evolution of pi (t) in time depends strictly 
on the concentration in the cell itself. 

The concentration of lAA in a cell depends on more factors than the pro- 
duction and the decay of auxin in the cell. Change of (t) is also determined 
by diffusion (passive transport) and active transport of auxin between the cells. 
The change over time of the concentration of lAA is modeled by the equation 



dtti (t) 



/^lAA^^w- E DKW-flj w) (2) 



dt 1 + K,^^a,{t) ^^^^ 

+ (ActiveTransportj_^j — Active Transport , 

where Pj^^ G K"*" is the lAA production coefficient, Kj^^ G K"*" is the coefficient 
which controls the saturation of lAA production, /ij^^ G M+ is the lAA decay 
constant and D G M+ is the lAA diffusion coefficient. The active transport 
depends on the presence of PINl denoted by pi and is modeled by the formula 

ActiveTransport,^, = T | W ^^^'""''l | °' (3) 

where T G M"*" is a polar lAA transport coefficient, h G is the exponentiation 
base which controls the extent to which the PINl protein distribution is affected 
by the neighboring cells and k,^ G is an lAA transport saturation coefficient. 
From equation ([2]) we know that the evolution of ai(t) depends only on itself, 
the first and the second neighbors of cell i. With the constructed graphs these 
neighbors can be determined easily since these are all the vertices in graph H 
with a distance to vertex i that is at most two. 

Equations ([I]), ([2| and ([3| describe the coupled model of Smith et al. that 
has been used to study the transport of hormones in the Arabidopsis shoot 
apex. In Smith et al| ( |2006[ ) they have studied a row of 50 equal sized cells with 



periodic boundary conditions. The time evolution starting from an initially flat 
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Figure 2: Simulation results of Smith et al. for a one dimensional row of 50 
cells with periodic boundary conditions and a two dimensional grid of irregular 
cells. The shades of green show a difference in the lAA concentration. The red 



shows the PINl concentration (figure reproduced from Smith et al (2006)) 



solution with a small amount of noise to break symmetry, showed the emergence 



of a pattern in the lAA concentrations (figure 2a I . Some cells have a very high 
concentration of the growth hormone. This results in a pattern with peaks that 
are equally spaced and become more prominent for an increasing lAA transport 
coefficient T. 

Similar patterns emerge in two dimensional models. There the cell shape is 
irregular. Again a pattern emerges with some cells having a high concentration 



of the growth hormone lAA (figure 2b ) 



2 The simulation problem 

2.1 The domain, its boundary conditions and the param- 
eters 

In this paper we analyze the solutions of equations ([T]) and ([2]), the coupled 
model of Smith et al., for a one dimensional file of equal sized square cells. This 
row of cells represents a part of the leaf from the left margin to the midvein. 
To provide the boundary conditions, two ghost cells are required at the end of 
the domain, since the model relates each cell with two cells at the left and the 



right (figure 3a) 



Th e graph and its weak dual for the one-dimensional domain are shown in 
fig 3b The length of each cell wall is equal to I. The set of vertices V of the 
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(a) 



-1 1 



n n+1n+2 



(b) 

Figure 3: (a) The one-dimensional model can be seen as a row of cells cut out 
of a leaf with equal sized square cells. In the first step we represent the leaf as 
a two dimensional squared grid of equal sized square cells. The second arrow 
indicates that we only consider the part from the left margin till the midvein. 
Here in each cell at the boundary, the direction of the auxin fluxes is indicated. 
In the third step a horizontal row of inner cells with, at each side two boundary 
cells is cut out of the domain. In the last step this domain is enrolled so it forms 
a one dimensional file of equal sized square cells. The two boundary cells at each 
side of the domain are the ghost cells, (b) The geometry of the leaf, used in 
the simulations, is a one-dimensional row of equal sized square cells. The weak 
dual is a row of connected cells. Each cell has a concentration of auxin that is 
denoted by ai{t) and a concentration of PINl that is denoted by Pi{t). The cells 
1 to n contain the unknowns, the cells —1, 0, n -I- 1 and n + 2 are ghost cells 
where the concentrations are set such that it produces in- and outgoing fluxes. 
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weak dual graph is split in n interior cells and four ghost cells, two on each side. 
The n interior points are labeled 1 to n. The ghost cells are cell —1 and at 
the left of the domain and cell n + 1 and n + 2 at the right. 

The concentration of the lAA hormone in the two ghost cells on each side 
are chosen to describe the influx at the boundary of the leaf and the efflux at 
the vein. The lAA concentration then changes linearly at the boundaries as if 
Neumann boundary conditions are applied. Until mentioned otherwise we will 
assume that the amount of influx is equal to the amount of efflux. This means 
that the boundary conditions become 

|a_i (t) ai (i) + 27 and ao (i) = ai (i) + 7, 
|a„+i (i) = a„ (i) - 7 and a^+i (t) ^ an (t) ~ 2-^ . 

The value of po (0 ^^^^ Pn+i (t) in the ghost cells is determined by equation ([I]) 
that couples it to the value of (t) in the ghost cell. Note that p^i (t) and 
^71+2 (t) do not appear in the problem since equation ([s]) does not require it. 
Together with an initial condition, the problem is transformed in an initial value 
problem that we can solve numerically with a time step method. Remark that 
these (homogeneous) Neumann boundary conditions are different from periodic 
boundary conditions. The concentrations in the cells on the left side of the 
domain can indeed be different from the concentrations in the cells on the right 
side of the domain. 

Equations p]) , (|2| and (|3| contain 11 parameters. A short description can 



be found in table 111 and further details can be found in Smith et al (2006). The 
values of these parameters must be real and positive. For the simulations in 
this paper we used three different parameter sets. Ml, M2 and M3. Parameter 



set M2 corresponds with the values used by Smith et al ( 2006 1 . Parameter set 
Ml and M3 contain the same values for the parameters as set M2 except for the 
lAA production coefficient. The value of this parameter is higher in parameter 
set Ml and lower in set M3. 

2.2 Time integration 

Similar to Smith et al, we can solve the initial value problem with numerical 
integration. Analysis of the eigenvalues of the Jacobian shows that they are 
mostly located along the negative real axis with some small complex conjugate 
outliers at the left of the imaginary axis. This suggests that the fourth order 



Runge-Kutta method (Hairer et al 2009) with time step At — 0.01 results in 



a stable method to integrate the equations. In figure |4] the time evolution is 
shown from i = to t = 10. The domain contains 20 cells plus 4 ghost cells 
where we assume zero Neumann boundary conditions. The parameter values of 
set Ml are used and the initial value for the concentration is 

Pi{t = 0)^5A and a,(i = 0) = 3.4, (5) 

where a small perturbation 0.2 sin ((5 {i + 2) tt) /24) for i = 1, . . . , 20 was added 
to it to break symmetry. 
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Table 1: Values of the parameters of equ ations (|l|)-([3| ) used in the simulations. 
The parameter values of M2 are found in Smith et al ( 2006 ) . 



Symbol 


Description 


Value 


Ml 


M2 


M3 


b 


Base for exponential PIN allocation 


3.000 


3.000 


3.000 


'''PIN 


PIN saturation coefficient 


1.000 


1.000 


1.000 




Transport saturation coefficient 


1.000 


1.000 


1.000 


'^lAA 


lAA saturation coefficient 


1.000 


1.000 


1.000 


PpiNo 


Base production of PIN 


0.000 


0.000 


0.000 


PpiN 


PIN production coefficient 


1.000 


1.000 


1.000 


MpiN 


PIN decay coefficient 


0.100 


0.100 


0.100 


MlAA 


lAA decay coefficient 


0.100 


0.100 


0.100 


PlAA 


lAA production coefficient 


1.500 


0.750 


0.500 


D 


lAA diffusion coefficient 


1.000 


1.000 


1.000 


T 


lAA transport coefficient 


3.500 


3.500 


3.500 




Figure 4: The time evolution of the solution for a row of 20 cells with zero 
Neumann boundary condition and parameter set Ml. The initial condition is 
given by equation ([5]) and we used RK4 for numerical integration. 



10 



Figure [4] shows the development of a pattern in the concentration of the growth 
hormone lAA. After a certain time the pattern arrives in a stable steady state. 
For the row of 20 ceUs a single peak with a high lAA concentration is formed. 



2.3 Steady state problem 

Rather than evolving the system in time, we can calculate the steady state 
solutions directly. We rewrite equations ([I]), ^ and ^ in order to obtain 
the steady state equations for this specific geometry and boundary conditions. 
Because the length of each cell wall is constant, I, it cancels from the equation 
since it appears both in the numerator and denominator of equation (|3| . The 
steady state problem becomes: 

„ PpiNn +PpiN^i • + • n 1 I 1 

= — MpinP* ^^^'^ I = 0,1, ... ,n,n + I 

1 + KpijjPi 

= MiaaOi - D (fli - ai_i) - D (fli - a,;+i) 



^, ft+ib°' \ ^( Kb°'+^ \ (6) 



b°' + b°*+^ ; 1 + K^af ^ Vb"'-' + b"*+^ ; 1 + K^af+i 

with i = 1, . . . . 



a_i = ai + 27 and cq = Oi + 7 
fln+i = an - 7 and a„+2 = - 27, 



where the indices and n + 1 in the first equation express the coupling of pi to 
ai in the first ghost cells. This system can be written as the system of equations 

F(U,A) = 0, (7) 

where F : M2n+m — ^ ^2n . A) F (U, A) with n the number of cells 
and m the number of parameters. U is a 2rt dimensional solution vector of the 
problem that contains both the p and a variables and A € denotes the set 
of parameters. 



3 The trivial solution 

In this section we search for a trivial solution of system (|6|, a solution we can 
calculate analytically that will be used as a starting point for the numerical 
continuation in section [5l 

If we assume that the solution is homogeneous and homogeneous boundary 
conditions are applied, i.e 7 = 0, then the values of p and a are the same for all 
cells so that 

Pi = pj and Ui = ttj Vz, J = — 1, . . . , n + 2. (8) 
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The system ([6| now reduces to 



= — — MpiN-Pi for « = 0, . . . , n + 1 

1 + KpiNl'i (9) 

0=— A^iaaO* tor z = 1, . . . ,n, 



and 

a„i = ao = ai and a„+i = a„+2 = a„. (10) 
Because and a,; are real positive numbers, the solution is given by 



-1 + ^/l + 4Kpj„ (ppj^^ + Ppj^fli) /Mpin 

2«:piN ' (11) 

lAA 



2k 



lAA 



with i = — l,...,ri + 2. This is the trivial solution of the system. 

From equation ([9| we know that for a certain parameter set, there is only one 
trivial homogeneous solution. By formula ( [Tl] ) it is easy to calculate this trivial 
solution for different parameter values. Figure [s] shows the concentration of lAA 
in one cell (cell number 6) versus the parameter Pj^A- Because the solution is 
homogeneous, the trivial solution curve would be the same for every cell and is 
independent of the number of cells. Figure [5] denotes also the trivial solution 
for parameter set Ml, M2 and M3 with a cross. 



3.1 Stability of the trivial solution 

We can calculate the trivial solution for every choice of the parameter, but, it 
will only be stable in a limited parameter range dependent on the number of 



cells. Although expression (11) for the trivial solution is simple, it is not easy 
to determine the stability of this solution. The Jacobian matrix of the coupled 
system ^ is not trivial and so are its eigenvalues. The stability can however 
be calculated by numerical means and in our simulations we approximate the 
Jacobian with central finite differences. The j-th column of J(U, A) is J(U, A)ej 
where is the unit vector with the j-th component 1 and the other components 
0. The column is then approximated as 

J(U,A)e, = ^(^ + -^'^)-^(^— (12) 

where e is taken of the order of 10~^. 

The stability of the trivial solution for a row of 20 cells is shown in figure 
[5j For smaller values of Pj^^ the eigenvalues of the trivial solutions lie in the 
left half-plane of the complex plane. Therefor these solutions are stable which 
is denoted with a full line in figure [5] For larger values of Pj^^ , at least one 
eigenvalue lies in the right half plane and so the trivial solution is unstable. 
This is indicated with a dotted line. 
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4 



< 
< 



Figure 5: The concentration of lAA in the trivial solution as a function of the 
lAA production coefficient Piaa- Each cell has the same concentration. However, 
for large Piaa this solution becomes unstable. Stability is calculated for a row 
of 20 cells. The other parameters are taken from table [TJ 

Also for other parameter values we can calculate the stability of the trivial 
solution. In each plot on figure [6] two parameter values are varied. The other 
parameter values are taken as in parameter set M2. On these plots we see for 
which values the trivial solution is stable (marked in gray) and unstable (not 



marked) for a row of 20 cells. For example figure 6a shows that for very small 
values of Pj^A almost every value of T gives a stable trivial solution. All other 
values of Pjaa give an unstable trivial solution if T is not too small. 

If we increase the number of cells, the shape of the stable region of the trivial 
solution remains approximately the same — it only gets slightly smaller. 



4 Methods 

4.1 Bifurcation analysis 

The study of the relation between the stability of a solution and the parameters 



of the corresponding dynamical system is known as bifurcation analysis ( Seydel 



19941. Such an analysis identifies the stable and unstable solutions and the 
bifurcation points that mark the transitions between them. This is biologically 
relevant since it will allow us to predict the patterns that emerge in the time 
evolution as the parameters of the model are changed. A bifurcation point is 
a solution (U,, A^) of system Q where the number of solutions changes when 
A passes A^. There are several types of bifurcation points, for example the 
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I- ,„ Q lO^^^r |_ 15 ,.- 

=L I W_ I \^ 

1 2 3 4 5 1 2 3 4 5 

P|AA P|AA D 

(a) (b) (c) 

Figure 6: Stable (marked in gray) and unstable region of the trivial solution for 
a row of 20 cells and different choices of the parameters Piaa (IAA production 
coefficient), D (IAA diffusion coefficient) and T (IAA transport coefficient). 
Other parameters are taken from parameter set M2. 



branch points, limit points and Hopf bifurcation points. A Hopf bifurcation is a 
transition where a periodic orbit appears and branch points and limit points are 
both bifurcation points among steady state solutions. A limit point, also called a 
turning point, is a point where, locally, no solutions exist on one side of the limit 
point and two solutions on the other side. A branch point on the other hand 
is a point where two or more branches with distinct tangents intersect. This 
analysis usually leads to a bifurcation diagram that highlights the connections 
between stable and unstable branches as the parameters change. It is useful to 
track all these solution branches that emerge, split or end in a bifurcation point. 
This can be done with the help of numerical continuation methods. 



4.2 Continuation methods 

The system of equations Q is a smooth map and we know that G Range {F). 
Following the implicit function theorem we know that for a regular point xq = 
(Uo, Ao) e M2n+m of F that satisfies F (xq) = 0, the solution set F'^ (0) can be 
locally parameterized about Xq with respect to some parameter s. This means 
that the system of equations F(U, A(s)) = defines an implicit curve U(A(s)) 



where A(s) : M — ^ M'" is any parametric curve in the M™ (Allgower and Georg 



19941. The idea of continuation methods is to find a curve c of approximate 
solutions U of the system in function of the parameter A(s). To construct such 
a curve of subsequent solution points x.; = (U^, A.^) — (U^, A(si)), continuation 
methods use a starting point xg = (Uq, Aq), a solution of system ([T]), along with 



an initial continuation direction (Krauskopf et al 20071. This starting point is 



typically a trivial solution. An important family of the continuation methods 
are the predictor-corrector methods such as pseudo-arclength continuation. The 
idea of the algorithm is to first predict a new solution point. In the corrector 
step, this predicted point is the start value for an iterative method that will 
approximate the solution to a given tolerance. For the pseudo-arclength, the 
predictor step uses the tangent vector to the curve at a solution point and a 
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given step size to predict a guess for the next solution point on the curve. The 
corrector step improves the guess with Newton iterations. 

Numerical continuation is available in AUTO (Doedel et al, 19971, LOCA 



part of Trihnos (Salinger et al 2005), PyDS (Clewley et al 2007) and others. 



These libraries can often also identify the bifurcations that occur along the 
continued curve and some of them, such as AUTO can automatically switch 
between branches at bifurcation points. 



5 Results 

In the next three examples we give, for a file of 20 cells, the numerical bifurcation 
analysis of problem ^ with respectively parameter sets Ml, M2 and M3. In 
the last example we also give the results of the bifurcation analysis of problem 
^ but now for a file of 100 cells instead of 20 cells and with parameter set Ml. 

In each example we impose homogeneous Neumann boundary conditions and 
the lAA transport coefficient T is the continuation parameter. Each time we 
display the bifurcation diagram that denotes the lAA transport coefficient T 
versus the lAA concentration in cell number 6. Remark that we can construct 
the bifurcation diagram with any measure on the y-axis. We have chosen T as 
the continuation parameter similar to the one-dimensional simulations of Smith 
et al. Also Jonsson et al. investigated the influence of the lAA transport 
coefficient T in their simplistic model by changing the ratio D/TP, with P the 
fixed value for PINl. 

We will find that the trivial solution loses its stability through either a branch 
point or a Hopf bifurcation. 

Example 1. The results of the bifurcation analysis for parameter set Ml 
are shown in figure [7j Figure [7a] shows the bifurcation diagram that depicts 
the concentration of the growth hormone auxin in cell number 6 versus the 
parameter T. The other plots in figure [7] show the steady state auxin patterns 
in all cells for the specific places indicated with labels in the bifurcation diagram. 
The trivial solution curve is the starting point of the continuation. It is the flat 
horizontal line in the bifurcation diagram. When the parameter T becomes 
larger than a critical value (T= 0.8983), the trivial solution loses its stability at 
a branch point. It was found by calculating for every solution (Ui, Ai) on the 
branch, the eigenvalues of the augmented Jacobian matrix deflned as 

J.ng ^ [JV \Jx] . (13) 

If the Jacobian Ju is singular and the rank of the augmented Jacobian is still 
smaller than 2n, then the solution point (Ui, A^) is a branch point. This means 
that there exist an eigenvalue fi (A.^) of the Jacobian which is equal to zero. 
Inserted into a graph, there is a path of an eigenvalue of the solution points 
corresponding to A close to A.^, that crosses the imaginary axis at the real axis 
when A = Aj. In the branch point on figure [7] there is an exchange of stability 



15 



4.0 



bifurcation diagram 




(a) Bifurcation Diagram 



2 10 



< 
< 



solution pattern 



5 10 15 

cell index 

(b) Solution 1 



solution pattern 




5 10 15 
cell index 

(c) Solution 2 



< 
< 



solution pattern 



5 10 15 
cell index 

(d) Solution 3 



.9 12 
§ 10 
" 8 



solution pattern 



5 10 15 

cell index 

(e) Solution 4 



solution pattern 



5 10 15 
cell Index 

(f) Solution 5 



solution pattern 




5 10 15 
cell index 

(g) Solution 6 



Figure 7: (a) The bifurcation diagram of example 1 for a row of 20 cells with 
continuation parameter T (lAA transport coefficient) versus the lAA concen- 
tration in cell number 6. Other parameters are taken from Ml. BP denotes a 
branch point and LP a limit point. The stars mark the places of the figures 
displayed below. 

(b)-(g) On these figures the lAA concentration in the whole domain is displayed 
corresponding with the stars marked on (a). 
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to another branch, also shown in the diagram. There are two stable parts on 
this other solution branch with patterns. When the lAA transport coefficient 
T is large, the stable solution pattern on this branch consist of one big peak 
(figure 7c I . The other stable part on this branch appears in a very limited range 
where T is smaller. For example solution 4 is such a stable pattern and it has 
two small variations (figure 7e I . The pattern in figure 7c is the same pattern 
that was obtained by numerical integration with the fourth order Runge Kutta 
in figure |4] as discussed in Section 2^ We thus found a connection between the 
trivial flat solution and the numerical solution with peaks. 



Example 2. The second example describes the results of the bifurcation anal- 
ysis for the model with parameter set M2 that Smith et al. used in their publi- 



cation (Smith et al 2006). The difference with the ffist parameter set Ml is a 
lower production coefficient of lAA, Pj^^. In the previous example, the stability 
of the trivial solution was lost in a branch point where branches with steady 
state solutions intersect. Now, we find that the stability is lost through a Hopf 
bifurcation where the equilibrium transitions into a periodic orbit. Looking at 
the eigenvalues of the Jacobian in this Hopf point, there is a pair of eigenvalues 
that satisfies 

^liX,)=±^|3. (14) 

If we draw a trajectory of the eigenvalues of solution points with A close to A^, 
we see that there are two complex conjugated eigenvalues different from zero 
that cross the imaginary axis when A = A^. 



Figure 8a shows the bifurcation diagram depicting again the concentration 
of lAA in cell 6 versus the continuation parameter T but now for parameter set 
M2. In this situation, the stability of the trivial solution is lost in a Hopf point at 
T = 3.3113. The branch that emerges from this Hopf point, shows the maximal 
and minimal lAA concentration over the orbit for each choice of parameter 
T. All the solutions on this branch are unstable and therefore we only have 
unstable periodic solutions. Further, also another steady state branch, different 
from the trivial solution branch, is displayed. This branch intersects with the 
trivial solution branch at a branch point (T = 5.4047 ). Around this branch 
point, all solutions are unstable. However, when we follow this new branch, we 
encounter another Hopf point where we now gain stability. The pattern of these 
stable solutions consist of one single big peak in the middle of the domain (see 
figure 8d| 



Example 3. For the third example parameter set M3 is used, which differs 
from sets Ml and M2 in the production coefficient of lAA. It is smaller than in 
set M2. 

The resulting bifurcation diagram is shown in figure [9j Just like in example 
2, the stability of the trivial solution is lost in a Hopf point (T — 22.7384 ). 
However, in contrast with this previous example, the periodic solution branch 
that intersects with the trivial solution branch in this point contains stable 
periodic solutions. Figure [TO] shows, in a three dimensional plot, the stable 
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Figure 8: (a) The bifurcation diagram of example 2 for a row of 20 cells with 
continuation parameter T (lAA transport coefficient) versus the lAA concen- 
tration in cell number 6. Other parameters are taken from M2. The stars mark 
the places of the figures displayed below. The dotted line through point HI 
shows the maximal and minimal value of the lAA concentration in cell number 
6 of the periodic solution for each choice of the parameter T. H and BP denote 
respectively a Hopf point and a branch point. 

(b)-(d) On these figures the lAA concentration in the whole domain is displayed 
corresponding with the stars marked on (a). 
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bifurcation diagram 




Figure 9: Bifurcation diagram of example 3 with parameter set M3. It de- 
picts the continuation parameter T (lAA transport coefficient) versus the lAA 
concentration in cell number 6. H denotes a Hopf point. The stable orbit for 
T= 23.5 is shown on figure 10 
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Figure 10: The time evolution for a row of 20 cells starting from the initial value 
in equation ( 15 ) with homogeneous Neumann boundary conditions, parameter 
set M3 but with lAA transport coefhcient T= 23.5. We used RK4 for numerical 
integration. The resulting solution is a periodic solution where the pattern 
changes from one peak concentration of auxin in the middle of the domain to 
a pattern with high concentrations at the boundaries. The periodic solution 
corresponds with the solution for T= 23.5 on the periodic solution branch in 
figure |9] 



periodic solution for lAA transport coefficient T = 23.5 found with RK4 starting 
from the initial value 

p,{0) = 1.79 and a,(0) = 3.76, (15) 

where a small perturbation 0.02 sin ((5 (i + 2) tt) /24) for i = 1, . . . , 20 was added. 
We see that the periodic solution changes in time from a pattern with one peak 
concentration of auxin in the middle of the domain to a pattern with two high 
auxin concentrations at the sides of the domain. In figure [TT] we plotted this tra- 
jectory in the (ag (t) , dag (t) /(ii)-plane starting from the initial value in equation 



Example 4. The previous three examples showed a part of the bifurcation 
diagrams corresponding with parameter sets Ml, M2 and M3, that differ in 
lAA production rate, for a one dimensional domain of 20 cells. In the next 
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Figure 11: The trajectory of the time evolution of figure 
(ae {t) , dae {t) /dt)-p\ane. We used a row of 20 cells starting from initial value 
( 15 ) with homogeneous Neumann boundary conditions, parameter set M3 but 
with lAA transport coefficient T= 23.5 and RK4 for numerical integration. 
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example we applied pseudo arc length continuation to a row of 100 cells and 
parameter values of set Ml. In this example, the stability of the trivial solution 



is again, as in example 1, lost at a branch point (T — 0.8504) (see figure 12a). 
The branch that crosses the trivial solution branch in this point is a bit more 
complicated. The branch contains 3 different stable parts. The stable part that 
contains solution point 2 consists of solutions with a pattern with 8 peaks (see 



figure 12c). Also the small stable area on the branch that contains solution 3 
consists of patterns with 8 peaks, but they are smaller due to the small value 
of T in this region. We see that the peaks become higher for an increasing lAA 



transport coefficient T (compare for example the patterns in figures 12c and 



12d or in figures 12f and 12g). The third stable part on this solution branch. 



appears also in a limited range of T. The patterns show only 7 peaks of high 



auxin concentration (see figure 12e). The patterns on the unstable part of the 



branch that contains solutions 5 and 6 also consist of 7 peaks (see figure 12f and 



12g 



The other solution branch in figure |12a| is found by time integration by 
starting with the steady state solution of figure |4] or the pattern in figure [7c] 
copied five times in a row which rapidly leads to a steady state that can be 
used as a starting point for the continuation. We see that this branch is not 
(directly) linked with the trivial solution branch and consists of one stable and 
one unstable part. On both parts the patterns consist of 9 high peaks of auxin 



concentration (see figures 12h and 12i 



Figures [Tj [Sj [9] and 12 for examples 1, 2, 3 and 4 show different bifurcation 



diagrams. In the examples 1 and 4, where parameter set Ml is used, the stability 
of the trivial solution is lost in a branch point. While in the examples 2 and 3 
it is lost in a Hopf point. We can calculate the type of these bifurcation points 
for every transition from stable to unstable for the three dimensional space of 
the parameters D, T and Pi^a- We found that in the one-dimensional case, for 
a row of cells, only these two different situations can occur: either the stability 
of the trivial solution is lost in a Hopf point, or the stability is lost in a branch 
point. Figure [T3| shows, for a row of 20 cells, parameter values from table [T] the 
parameters D, T and /Ojaa) the corresponding type of bifurcation. For example 
figure |13a| shows that for small values of the production coefficient of lAA the 
stability of the trivial solution will be lost in a Hopf point. Indeed, in example 
2 and 3 we found similar results. 



6 Conclusion and discussion 



In this paper we have explored the model proposed by Smith et al (2006) for 



the transport of growth hormones in a one-dimensional row of cells. This model 
describes the evolution of the concentration of PINl and lAA in each of the 
cells by a coupled set of non-linear ordinary differential equations. The change in 
concentration of lAA in each cell depends only on PINl and lAA concentrations 
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Figure 12: (a) Bifurcation diagram of example 4 for a row of 100 cells and param- 
eter set Ml. The diagram shows the continuation parameter T (I A A transport 
coefficient) versus the lAA concentration in cell number 6. The unconnected 
branch is found by time integrating the steady state solution of figure [7c] copied 
five times in a row until a steady state is found. The resulting steady state is 
then used as a starting point for the continuation. The stars mark the places of 
the figures displayed below, (b)-(i) On these figures the lAA concentration in 
the whole domain is displayed corresponding with the stars marked on (a). 
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(c) p,^^ = 0.750 

Figure 13: Map of types of bifurcation points for a row of 20 cells and different 
choices of the parameters Pj^^ (lAA production coefficient), D (lAA diffusion 
coefficient), T (lAA transport coefficient). Other parameters are taken from 
table [1] 
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in that cell, its nearest and next-nearest neighbors and this leads to a sparsely 
coupled system. 

We have analyzed the steady state solutions of the system as a function of 
three of the 11 parameters in the problem. We have varied the lAA transport 
coefficient T, the diffusion coefficient D and the production rate Piaa- We have 
used numerical continuation to generate these solutions starting from a trivial 
solution of the system. 

The trivial solution is identified as an analytical solution where the concen- 
tration in all the cells is the same. This solution is stable for some region of the 
parameter space. However, changing the parameters, for example, increasing 
the lAA transport coefficient T, destroys its stability. In contrast to the uncou- 
pled system studied by Jonsson et al (2006), the eigenvalues of the Jacobian of 
this coupled system can not be easily analyzed analytically. The Jacobian now 
has a blocked sparse structure and its eigenvalues are studied numerically. 

In the exploration of the solutions, we have identified two generic bifurcation 
scenarios through which the trivial solution loses its stability. These scenarios 
reappear for various choices of the parameters. In the first scenario, a stable 
solution can lose its stability through a branch point, where it becomes a pattern 
with regular spaced peaks of high auxin concentration. This solution was already 
found by Smith and collaborators by direct numerical simulation. The spacing 
and the height of the peaks in the pattern depends on the other parameters of 
the system. 

However, we have found in the coupled system that the trivial state can 
also lose it stability through a Hopf bifurcation, where the Jacobian has two 
complex conjugate eigenvalues that become purely imaginary. For a limited 
parameter range this leads to stable periodic solutions, where the concentration 
in the cell changes periodically over time. However, for another range of pa- 
rameters these periodic orbits are unstable and the trivial solution then loses its 
stability through an unstable orbit. The steady state solution then falls, beyond 
a parameter threshold, back to a pattern of regularly spaced peaks (see figure 
[S]). These Hopf bifurcation and the periodic solutions are not present in the 
model studied by Jonsson where the PINl concentration is kept constant and 
all eigenvalues of the Jacobian are real. 

These are the only two bifurcation scenarios that we have found at the 
stability boundary of the trivial solution for various choices of the 11 parameters 
in the model and for an increasing number of cells in the row. 

Although the paper studies the steady state solutions of a rather academic 
model with a row of equal sized square cells, the authors believe it is a valuable 
contribution to our understanding of pattern formation by auxin accumulation 
since it builds the foundation for a rigorous bifurcation analysis of the steady 
state patterns in a two dimensional array of cells. We expect to find there a 
similar trivial solution that will lose again its stability as the parameters change 
and turn into regular patterns of high concentration peaks and time periodic 
solutions. 

In this paper we have studied bifurcations with 7 = 0, so called homogeneous 
Neumann boundary conditions. Similar studies can be done for inhomogeneous 
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Neumann boundary conditions, however, then there is no analytically solvable 
trivial solution. We have explored continuation in 7 and found, amongst others, 
s-shaped bifurcation diagrams with double limit points. This leads to a hys- 
teresis effect in the boundary conditions. This will be further explored in the 
future. 

There are still many uncertainties in the current generation of models. Es- 
pecially the large number of parameters and the uncertainty in their values is a 
reason for concern. By focusing on the qualitative properties of the transitions 
that appear in the models rather than on the states for particular choices of 
the parameters, we hope to understand more about the possible patterns that 
appear in real systems. It is valuable to calculate similar bifurcation diagrams 
for all the proposed models for auxin transport using the numerical continuation 
methods. This will allow the comparison of models across a range of parameters 
and check if they exhibit qualitatively the same transition if parameters change. 

In real plants it is impossible to tune a parameter such as the transport 
coefficient T in a continuous way as is done in these calculations. It can only 
be changed in discrete steps in a plant by the introduction of, for example, mu- 
tations that compromise or enhance the auxin production or transport. Com- 
paring such experiments with the model outcome will make it possible to refine 
the models to give a more realistic description of the biological system. 

The bifurcation analysis also yields interesting new insights into the potential 
behaviors of the biological system. It is interesting to note that low values of the 
lAA transport coefficient lead to flat distributions of the auxin concentration, 
whereas high concentrations are required to establish sharp accumulation peaks. 
Indeed, experimental inhibition of auxin transport with NPA abolishes the nor- 
mal narrow accumulation peaks and results in a much flatter auxin distribution 
pattern ( Scarpella et al 2006 ) . 

One especially interesting behavior is the oscillation obtained with a specific 
set of parameter values (Fig 11). This behavior has to our knowledge never 
been observed in the context of leaf development. However, in the root basal 
meristem, oscillating auxin concentrations have been observed and related to 



the regular induction of laterals along the growing axis of the root (De Smet 



et al 20071. 



One other characteristic unveiled by the bifurcation analysis is that across 
the region where a pattern of auxin accumulation peaks are generated, these 
peaks occur at very stable distances. This would imply that in the case of vas- 
cular patterning, the initial distance between veins is relatively stable, implying 



that observed differences in vascular density in mature leaves (Dhondt et al 



2011 1 largely result from differences in subsequent development. This is an ex- 



ample of an new hypothesis generated by modeling a biological process that 
can be experimentally validated and underlines the importance of systematic 
exploration of biologically relevant parameter variations. 

It is important to repeat that we have kept the plant geometry fixed in the 
current model. It is an open question how the calculations can be extended to 
include cells to undergo growth and division. 
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